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Abstract 

We describe a method to calculate the electrical force acting on a sphere in a suspension of 
dielectric spheres in a host with a different dielectric constant, under the assumption that a spatially 
uniform electric field is applied. The method uses a spectral representation for the total electrostatic 
energy of the composite. The force is expressed as a certain gradient of this energy, which can be 
expressed in a closed analytic form rather than evaluated as a numerical derivative. The method is 
applicable even when both the spheres and the host have frequency-dependent dielectric functions 
and nonzero conductivities, provided the system is in the quasistatic regime. In principle, it 
includes all multipolar contributions to the force, and it can be used to calculate multi-body as 
well as pairwise forces. We also present several numerical examples, including host fluids with 
finite conductivities. The force between spheres approaches the dipole-dipole limit, as expected, at 
large separations, but departs drastically from that limit when the spheres are nearly in contact. 
The force may also change sign as a function of frequency when the host is a slightly conducting 
fluid. 
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I. INTRODUCTION 



An electrorheological (ER) fluid is a material whose viscosity changes substantially with 
the application of an electric field 1]. Generally, such fluids are suspensions of spherical 
inclusions of dielectric constant ti in a host fluid of a different dielectric constant eh- The 
viscosity is believed to change because the spheres acquire electric moments (dipole and 
higher) when an electric field is applied, then move under the influence of the electrical forces 
between these induced moments. These forces typically cause the spheres to line up in long 
chains parallel to the applied field, thereby increasing the viscosity of the suspension. The 
viscosity relaxes to its usual value when the field is turned off, and the chain-like structure 
disappears. 

ER fluids have potential applications as variable viscosity fluids in automobile devices^], 
vibration control(3|, and elsewhere. Furthermore, their operating principle is also relevant 
to other materials, such as magnetorheological (MR) fluids 4]. These are suspensions of 
magnetically permeable spheres in a fluid of different permeability, whose viscosity can be 
controlled by an applied magnetic field. 

To obtain a quantitative theory of ER (and MR) fluids, one needs to understand the 
electric-field-induced force among the spheres. At low sphere concentrations and large inter- 
sphere separations, this force is just that between two interacting electric dipoles whose 
magnitude is that of a single sphere in an external electric field. But at smaller separations, 
the force deviates from the dipole-dipole form. Besides this electrostatic interaction between 
the spheres, there are other forces acting on the spheres, including a viscous frictional force 
from the host fluid, and a hard-sphere force when the two dielectric spheres come in contact. 
In the present paper, we will be concerned only with the electrostatic force. 

A number of existing theoriesgo beyond the dipole-dipole approximation in calculating 
electrostatic forces in ER ffuidsjjy, Q, 0, IQ, U, llE 13 5 and several experiments have 



)een carried out which are relevant to forces in the non-dipole regime (see, e. g., Refs. 
3 13)- Klingenberg et al^\ have incorporated both multipole and multi-body effects into 

I I 

the sphere-sphere interactions, using a perturbation analysis. Chen et a/. [6] have described 
a multipole expansion for the forces acting on one sphere in a chain of spheres in a fluid 
of different dielectric constant, and find a strong departure from the dipolar limit when the 
particles are closer than about one diameter. Davis p] has calculated the electrostatic forces 



3 



between dielectric spheres in a host fluid directly, using a finite element approach to solve 
Laplace's equation for a chain of particles in a host dielectric. In a more recent work^, 
he has used an integral equation approach to calculate the interparticle forces in ER fiuids, 
including effects due to time-dependent application of an external field^and nonlinear fiuid 
conductivity. A finite-element approach has also been used by Tao et al. j9| to solve Laplace's 
equation and obtain the electrostatic interactions between particles in a chain of dielectric 
spheres in a host fiuid; they found, as in Ref. [3], that the dipole-dipole approximation is 
reasonably accurate for large separations or moderate dielectric mismatches, but fails in 
closely spaced particles and large mismatches. Clercx and Bossis|3| have gone beyond the 
approximation of dipolar interactions to include multipolar and many-body interactions, 
expressed in terms of the induced multipole moments on each sphere; they also obtain an 
expression for the forces in terms of these induced multipole moments. 

As discussed further below, the electrical force acting on a sphere in an electrorheological 
fiuid is basically the gradient of the total electrostatic energy of that fiuid with respect to 
the position of the sphere. This total electrostatic energy can, in turn, be expressed in 
terms of the effective dielectric tensor of the suspension, a quantity which has been studied 
since the time of Maxwell. Indeed, numerous authors have calculated this tensor in a wide 
variety of geometries, going well beyond the regime of purely dipolar interactions. For 
example, Jeffrey 3] has calculated the total energy of two spheres in a suspension as a 
function of their separation and the dielectric mismatch. From this total energy, the force 
can be obtained numerically as the derivative of this energy with respect to separation. 
Recently, the pairwise forces between spheres of different sizes have been calculated using 
the so-called dipole-induced-dipole approximation, and even approximately including the 
effects of other spheres 3l- Once again, the forces were obtained explicitly by numerically 
differentiating the total electrostatic energy with respect to particle coordinates. McPhedran 
and McKenziej3], and Suen et a/.[3l, and many others, have calculated the total energy 
of spheres arranged in a periodic structure. In principle, forces could also be extracted 
from this calculation by taking numerical derivatives, provided that the distortions of the 
structure leave it periodic. The energy of a non-periodic suspension of many spheres has 
been studied by Gerardy and Ausloosj2^ and by Fu et al. 21 1, in both cases including 



large numbers of multipoles. Once again, forces can be extracted, in principle, from these 
calculations by taking numerical derivatives of the computed total energies with respect to 
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sphere coordinates. 

Several authors have included the effects of finite conductivity on forces in electrorhe- 
ological fluids, and have also considered how such forces depend on frequency. Davis j^] 
has analyzed polarization forces and related effects of conductivity in ER fluids. Tang et 
g/.jiol have calculated the attractive force between spherical dielectric particles in a 
conducting film. Khusid and Acrivos^^ have considered electric-field-induced aggregation 
in ER fluids, including interfacial polarization of the particles, the conductivities of both 
the particles and the host fluid, and dynamics arising from dielectric relaxation. Claro and 
Rojasj3] have calculated the frequency-dependent interaction energy of polarizable particles 
in the presence of an applied laser field within the dipole approximation; they considered 
primarily optical frequencies rather than the low frequencies more characteristic of ER fluids. 
Ma et al. 241 have considered several frequency-dependent properties of ER systems, starting 
.o,n a weuinow. spec.a. .ep.ese.ta^o JflfHQfi6 fo. the die.ectnc f..ti„n 



of a two-component composite medium. Finally, Huang |31l| has carried out a calculation of 
the force acting in electrorheological solids under the application of a non-uniform electric 
field, and considering both finite frequency and finite conductivity effects. 

A common feature of most of the above approaches is that they involve first calculating 
the total electrostatic energy of the suspensions, then obtaining the forces by numerically dif- 
ferentiating this energy with respect to a particle coordinate. This numerical differentiation 
is cumbersome and can be inaccurate. Ref. [3| does give an expression for the force, but 
in terms of implicitly defined multipoles. In this paper, by contrast, we describe a method 
for calculating these forces explicitly, without numerical differentiation. This approach is 
computationally much more accurate than numerically differentiating the energy. While our 
new method may appear to be merely a computational advance, its additional accuracy and 
flexibility should make it widely applicable. 

Specifically, our approach allows one to calculate the electric-field-induced force between 
two dielectric spheres in a host of a different dielectric constant, at any separation. It is 
applicable, in principle, to spheres of unequal sizes, to particles of shape other than spheres, 
to suspensions in which either the particle or the host or both have nonzero conductivities, 
and to systems whose constituents have frequency- dependent complex dielectric functions. 
It can also be used to calculate the electrostatic force on one particle which is part of a 
many-particle system, and thus is not limited to two-body interaction. It should thus be 
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useful in quite general circumstances including, in particular, non-dilute suspensions. 

Our approach starts, as do previous calculations, with a method for calculating the total 
electrostatic energy of a suspension of (two or more) spheres in a host material of different 
dielectric constant. We choose to express this total energy in terms of a certainpole spectrum 



arising from the quasistatic resonances of the multi-sphere svstem|2c 



26 



This representation has previously been used to calculate the f req uency- de p en dent shear 



modulus, static yield stress, and structures of certain ER systems 



33, M . The force 



on a given sphere in a multi-sphere system involves a gradient of this energy with respect 



to the positiori o 
previous work 



that sphere. But rather than evaluating this derivative numerically, as in 
24| . we express this derivative in closed analytical form in terms of the pole 
spectrum and certain matrix elements involving the resonances. This expression is readily 
evaluated simply by diagonalizing a certain matrix, all of whose components are readily 
computed. 

Our approach has formal similarities to the well-known Hellmann-Feynman expression 
for forces in quantum-mechanical systemsjs^. In the quantum-mechanical case, the force is 
expressed as the negative gradient of system energy with respect to an ionic position. This 
energy is the expectation value of the Hamiltonian in the ground state. According to the 
Hellmann-Feynman theorem, the gradient operator can be moved inside the matrix element, 
thereby eliminating the need to take numerical derivatives. The Hellmann-Feynman force 
expression is the basis for many highly successful molecular dynamics studies in quantum 
systems (see, e. g., Ref. j^)- present classical case, the total energy can also be 

expressed as a certain matrix element of an operator, and thus, just as in the quantum 
problem, the force is the gradient of that expectation value. In this paper, we shall show 
that, again as in the quantum case, the gradient operator can be moved to within the matrix 
element, and the need to take a numerical derivative is eliminated. This simplification allows, 
in principle, the calculation of forces in very complicated geometries, even though, in the 
present paper we shall give only relatively simple numerical illustrations involving forces 
between two spherical particles. 

The remainder of this paper is organized as follows. In Section |Hl we present the for- 
malism necessary to calculate the forces in a system of two or more dielectric spheres in a 
host medium, without taking a numerical derivative. In Section ITTTl we give several numer- 



ical examples of these forces, at both zero and finite frequencies, for a two-sphere system. 
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Section HVl presents a concluding discussion and suggestions for future work. 



II. FORMALISM 

Let us assume that we have a composite consisting of spherical inclusions of isotropic 
dielectric constant ei in a host of isotropic dielectric constant eh, both of which may be 
complex and frequency- dependent. We will assume that a spatially uniform electric field 
RefEoe"*"^*] is applied in an arbitrary direction (we take Eg real). We also assume that the 
system is in the "quasistatic regime." In this regime, the product A;^ <^ 1, where k is the 
wave vector and ^ is a characteristic length scale describing the spatial variation of e(x, cj). 
Under these conditions, the local electric field E(x, cj) = — V$, where $ is the electrostatic 
potential. Finally, we assume that the R*'^ spherical inclusion is centered at R, and has 
radius or. The approach which we use automatically includes all local field effects. 

Since our force expressions differ slightly at zero and finite frequencies, we will first present 
the formalism at = 0, and then generalize the results to finite uj. 

A. Zero Frequency 

If the position of the spheres is fixed, the total electrostatic energy may be written in the 
form 

w=^j:j:^o,,eo,eo,, (1) 
^ ^ i=i j=i 

where V is the system volume, ee-^ij is a component of the macroscopic effective dielectric 
tensor, and i?o,i is a component of the applied electric field. Eq. ((H) is, in fact, a possible 
definition of ee;ij(co') Q|. To produce this applied field, we require that $(x) = — Eq-x at the 
boundary S of the system, which is assumed to be a closed surface enclosing V. In writing 
eq. (^, we allow for the possibility that the spheres in the composite are arranged in such 
a way that the composite is anisotropic even though its components are not. 

For an iso trop ic composite, ep mav be written in terms of a certain pole spectrum of the 



ijroptc composite, epU 



composite as 
where 



Sa is a pole, and Ba is the corresponding residue. The poles Sa are confined to the interval 
< < 1. For an anisotropic composite, this form may be generalized to 

eh 

where 6ij is a Kronecker delta function and Ba-ij is a matrix of residues. This form is 
general, applicable to any two-component composite material which is made up of isotropic 
constituents, but is not necessarily isotropic macroscopically. As in the isotropic case, the 
poles are confined to the interval < Sq, < 1. 

The poles Sa are the eigenvalues of a certain Hermitian operator F, and the residues Ba 
are determined by the eigenvectors of that operator. F is defined in terms of its operation 
on an arbitrary function 0(r) by the relation 



F0(r) = / d^r'Vl-^—] ■ V'0(r') 
J vtot y I r — r \ J 



(5) 



where the integration runs over the total volume vtot of all the inclusions of ei. As in Ref. 
37[, we introduce a "bra-ket" notation for two potential functions to denote their inner 
product, 

{(j)\ij)=f dVV0*(r) ■ V^(r). (6) 

"'ftot 

Physically, the eigenvalues correspond to the frequencies of the natural electrostatic modes 
of the composites, at which charge can oscillate without any applied field, and the corre- 
sponding eigenvectors describe the electric fields of those modes. 

It is convenient to express F in terms of its matrix elements between the normalized 
eigenstates of isolated spheres. In this basis, and using eqs. © and ©, it is found that F 
has the following matrix elements: 

FR^m;R'£'m' = {i'mml^ IpK'e'm') = Si 5i^£' Sm,m' Sr^B.' + <5R£m;R'£'m' (1 " (^R.R')' (''') 

where S£ and Q-Rem;ii'e'm' will be given further below. Inside the sphere centered at R, 
'j/'Rfmli") is equal to an eigenfunction or resonance state of that isolated sphere, while outside 
that sphere ipnimi''^) = 0. The angular dependence of ipnimi''^) is given by the spherical har- 
monic YimiO, 0), which has an order im multipole moment of electric polarization. However, 
the eigenvalue of this state depends only on £: 

£ 

S-Rim = Si= 2£j^l- 
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The quantity Qn£m;'R'e'm' represents the matrix element of T between two states of two 
different, non-overlapping spheres (i. e., |R' — R| > + or/) and is given by 

^+1/2 £'+1/2 / U/2 



X 



R' - R|^+^'+i V (2 £ + 1) (2 f + 1) ^ 
{i + i' + m- m')\ 



[{i + m)\{e - m)!(f + m')!(f - m')\]^^^ 
X e*^i^'-i^("'-"^)p^7;V'"(cos^R/_R), (9) 

where ^r'_r and 0r'-r are polar and azimuthal angles of the vector R' — R, and the 
functions P^T^^™ associated Legendre polynomials. 

If we denote the eigenfunctions of T by ipa{r), then Sa and '?/'a(r) satisfy the eigenvalue 
equation 

TtP^{r) = sMr). (10) 

Since F is a Hermitian operator, the eigenvalues Sa are real, and the corresponding eigen- 
functions are orthogonal and can be chosen to be orthonormal. Again, it is convenient to 
represent them using a bra-ket notation. In this notation, the eigenfunctions are denoted 
la) and the orthonormality condition is 

{a\a') = Sa,a' ■ (11) 

The eigenvalues Sa are the poles of eq. or eq. (jH). 

The corresponding residues B^^ij may be expressed in the same bra-ket notation as 

B^.,, = ^-^{t\a){a\j)^M:Mr, (12) 

The matrix element = {i\a) is basically the component of the electric dipole moment of 
the eigenfunction \a) in the i^^ Cartesian direction. 

It is convenient to expand both the eigenfunctions |a) and the states \i) {i = x, y, z), in 
terms of the single-sphere eigenfunctions '?/'R£m(r) mentioned above. In bra-ket notation, 

l«) = E A^mm^)- (13) 
Rem 

The expansion coefficients satisfy the normalization condition J^Rem l^mmP — ^5 where the 
indices i = 1,2,... and m = —i, respectively. Similarly, the states \i) may 

be expanded as 

) = E M^J^^^)^ (14) 



Rim 



where i = x, y, z. If the z axis is chosen as the polar axis for the spherical harmonics, then 



the M-^g^ take the formjgel 



V ^ "^tot / 



1 /2 

MA,m = (—) S^,oki- (15) 



Thus, the matrix elements are given explicitly by 

= E(|^f'(^Rll + ^Rl-l), 



R 



m: = E {yf' ^Rio- (16) 

In other words, the residues of the a*^ eigenfunction are basically the square of the electric 
dipole moment of that mode in the x, y, or z direction. 

Combining these results, we can re-express the matrix elements ((3]) of the dielectric tensor 
in bra-ket notation first as 

r ee;ij _ vtot ^ {i\a){a\j) . . 

where the explicit forms of and {ot\j) are given by eqs. 

We now use the above formalism to obtain an expression for the force on a dielectric 
sphere centered at R in a suspension consisting of an arbitrary assembly of spheres. First, 
we rewrite eq. (fTTj) as 

5.,-^ = ^(^|G(.)|j), (18) 



where 



^ ^ MM = (,j _ rr\ (19) 



is a Green's function for this problem, I is the identity matrix, and the matrix elements of 
r are given by eq. ((Tj). If the applied electric field is Eq, the total energy takes the form 

= ^ eh Eo ■ Eo - ^ E Eo, (z|G(.) |j) E,,. (20) 
o vr o TT . . 
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We now write the k component of the force on the sphere at R as 

---(HI 

Here Rk denotes the A;*'^ component of R, and the subscript $ denotes that the derivative 
is taken with the potential fixed on the boundaries. The positive sign, though seemingly 
counterintuitive, is actually correct here because the system is held at fixed potential on the 
boundaries 13]. Using eq. (plj) . the derivative in eq. (PT|) can be expressed as 

'dW^ _ 



:I-'^?^"'-^-<''feH'^>^ '''' 



\dRk 

The derivative can be brought inside the bras and kets because these bras and kets do not 
depend on Rk. 

The derivative of G{s) appearing in eq. ()22|) can be evaluated straightforwardly. Let us 
assume that the operator F depends on some scalar parameter A (e. g., Rk)- Then, if we 
introduce the operator Ux = dT/dX, we can calculate the partial derivative dG{s, X)/dX as 
follows: 

dG{s,X) 

lim 



dX dx^o 
= lim 

dX^O 



{sl - r(A) - Ux dX}-^ - {si - r{x)}-^\ I dX 
{I - {si - r(A))-i Ux dx}-\si - r(A))-^ - {si - r(A)}-i] / dX 
= [si-r{X)]''Ux[si-r{X)]-' 

= G{s,X)UxG{s,X). (23) 

We can now use the above identity to calculate the force as given in eqs. fET|) and ((221) ■ 
The result is 

Fr, = -^J:J:Eo,E,,{z\G{s,X)U^,G{s,X)\j), (24) 

* j 

where we have introduced 

Using the representation (fTIHl for G{s, A) [and taking the eigenvalue Sa and the eigenstate 
la) to refer to the operator T{Rk)], we can rewrite eq. as 

Stt ''^y is-s^)is-Sf,) ^ ^ 

Eq. ()26|) is our central formal result. 
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As noted earlier, eq. ()22j) bears a resemblance to the Hellmann-Feynman theorem in 

IT" 

quantum mechanics |35]: in both cases, the derivative of an operator with respect to a 
parameter appears inside a matrix element. But there is a significant diff'erence between 
the two. In the Hellmann-Feynman case, the ket which plays the role of \i) is an eigenstate 
of an operator, which is the actual Hamiltonian of the system. Although the ket in that 
case depends on A, the derivative can still be moved inside the bra and ket because the 
eigenstates are orthonormalized. Here, by contrast, the states \i) and |j) are not eigenstates 
of the operator F, but they do not depend on A; so the derivative can still be moved inside 
the matrix element. This simplification allows forces to be computed without carrying out 
numerical derivatives. 

Eq. (j26|) may appear to be rather difficult to apply in practice. But in fact it is compu- 
tationally quite tractable. Basically, there are two matrices which are needed as inputs: G 
and Uji^. = dr/dRk- G is diagonal in the same basis as F. All the matrix elements of F are 
explicitly known in the Him basis [cf. eq. ((Zj)]. Likewise, the matrix elements of dT/dRk 
can be obtained from F purely by elementary calculus. Thus, in order to compute the force 
component Fr^, one first finds the eigenvalues and eigenfunctions of F (and hence of G), 
then computes the matrix GUr^^G in the basis in which F is diagonal, and finally the ma- 
trix elements {i\G U iif,G\j) , from which the force can be computed for any direction of the 
applied field Eq. Since the diagonalization can be done with standard computer packages, 
the whole procedure is well defined and straightforward. Furthermore, once F has been 
diagonalized, the same basis can be used to compute the forces for any value of the variable 
s = {l-e-Je^)-\ 

To illustrate how this formalism can actually be used to compute the force explicitly, 
we will consider just a suspension of two spheres, the two spheres being located at (0, 0, 0) 
and (0, 0, Rq). The total energy is given by eq. ((H). We consider two configurations for the 
electric field: Eq = (0,0, i?o) (we call this the "parallel configuration") and Eq = (-E'o,0,0) 
("perpendicular configuration"). In both cases, the component of the force on the sphere at 
Rq along the axis joining the two spheres can be calculated using eq. ()2(i|l . 

To compute the force explicitly in this example, we have to consider how the operator F 
changes with the separation Rq of the spheres, so that we can compute the matrix elements 
of Ur^. According to eqs. ^ and 0, the diagonal matrix elements of F are independent of 
Rq, while each of the off-diagonal matrix elements, according to eq. 0, is proportional to an 
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integer power of 1/|R' — R| = 1/Rq. Hence, Ur^, = dV j dRo is easily calculated in a closed 
form. For the case of two spheres, it is straightforward to calculate this derivative. The 
eigenstates \a), as well already known if the original eigenvalue problem involving 

r(i?o) has been solved. The kets \i) are given by eq. (jl5j) . Therefore, it is straightforward 
to calculate the quantity and hence the force, using eq. ^2^i. 

B. Finite Frequencies 

The results of the previous subsection are readily generalized to finite frequencies. In 
this case, the total electrostatic energy will be a sinusoidally varying function of time. The 
quantities of experimental interest will be the time-averaged electrostatic energy W^v and 
time-averaged forces. W^v is given by the generalization of eq. IQ, with an extra factor of 
1/2 to take into account time- averaging, namely 



V 

16^ 



Re 



(^e;ijiuj) Eo^iEoj 



i=l j=l 



(27) 



Here the applied field is assumed to be Eq cos(ciJ)f:) = Re[Eo e"*"^*], ee-ij{uj) is a component 
of the complex frequency-dependent macroscopic effective dielectric tensor, and Eq is a real 
vector. All the remaining equations in Sec. ITTK continue to be valid up to eq. ()21|). which is 
replaced by 



+ 



\ 9Rk 

where $ is held fixed. The generalization of eq. (^^1) is 

3 3 



av,_Rfe 



-Re 



16 vr 



{z\a){a\UR,\P){P\j) 



i=i j=i 

Expression can be evaluated just as at u 
frequency can also be computed explicitly. 



(28) 



(29) 



0, and thus the time-averaged force at finite 



III. NUMERICAL RESULTS 



We have applied the above formalism to two spheres of dielectric constant e; in a host 
of dielectric constant eh. In most cases, we assume that the spheres have the same radius. 
We choose a coordinate system such that the two spheres are located at the origin and at 
R = -Rz, and we consider two configurations for the applied electric field, Eq = Eqz and 
Eq = EqU, as shown in Fig. [TJ 
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Once the elements of the T and Ur matrices are known, the calculation of the interparticle 
force reduces to an eigenvalue problem. To carry out the various required matrix and vector 
operations, we used GNU Scientific Library (GSL) routines^] and C++ complex class 
library. In the parallel configuration, we calculated all the elements in the F and f//j matrices 
up to £max = 80; it is easy to include such a large cutoff because only m = needs to be 
considered for this geometry, the polar and azimuthal angles of R equaling zero. Despite the 
large cutoff, most of the contributions to these matrices came from i < 10. Based on this 
information, we set £max = 10 for the F and Uji matrices in the perpendicular geometry. Even 
with this cutoff, the matrices involved in this calculation are large since m can be nonzero 
in the perpendicular case: the dimension of the matrix for £ < 10 is 2X]l=i(2^ + 1) = 240. 

The F matrix for both cases consists of four square blocks. The two diagonal square blocks 
have diagonal elements = i/{2i + 1) with all off-diagonal elements vanishing. The other 
two (off-diagonal) square blocks have elements Qoim^iu'm'- For the Ur matrix, the diagonal 
square blocks have all zero elements, and the elements of the two off-diagonal blocks are 
equal to dQo£m;m'm' / dR. Once we have calculated all the eigenvalues and eigenvectors of 
the F matrix, we can compute the Ma, Ba, and hence the force on the sphere from {a\UFi\(3), 
using eq. or (OHl). 

As a first example, we have considered ej = 10^, Ch = 1. The choice for ei approximates 
the value ei = oo corresponding to two metallic spheres at zero frequency in an insulating 
host with unit dielectric constant. In Fig. |21 we show the magnitude of the calculated radial 
component of the force acting on the sphere at R, as a function of the sphere separation, for 
both parallel and perpendicular configurations. Although not apparent from the plot, this 
component of the force is attractive (i. e., negative) in the parallel configuration, repulsive 
in the perpendicular configuration. We have arbitrarily chosen sphere radii of a = 3.15 mm 
and a field strength of Eo = 25.2 V/ mm as in recent experiments carried out in Ref. Q| 
(for different materials). However, the forces are easily scaled with both field strength and 
sphere radii: for fixed €[ and eh the appropriate scaling relation is 
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a%V±,||(ei,eh,i?/a), (30) 



where f± and /y are functions of ej, eh and the ratio R/a. 

It is of interest to compare these plots with the same forces as calculated in the dipole- 
dipole approximation. For two parallel dipoles pi and p2, located at the origin and at 
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R = i?z, the z component of the force acting on the sphere at R has the well-known form 

-^^P'"(i?) = 3^fl-3(pi-zfl, (31) 
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where pi and p2 are the magnitudes of the two dipole moments, and pi is a unit vector 
parallel to pi (or P2). For the present case, if the spheres are well separated and have equal 
radii, the dipole moments can be calculated as if each is an isolated sphere in a uniform 
external electric field Eq: 

Pi = P2 = a'Eo^^. (32) 

For the cases in which the unit vector Eg is perpendicular and parallel to z, the radial 
component of the force reduces to 



12 — ~2 12 — OU l^Q 



_ei + 2j 

These values of -^1^2'''''' and F^g''" shown in Fig. |21 agree very well with those calculated from 
eq. (jH^ (taking ei = 00) at large separation {R ^ a) but depart strongly at small separation 
{6 = R — 2a <^ a). Just as in the exact calculation, the radial component of the force in 
the dipole-dipole limit is repulsive in the perpendicular case and attractive in the parallel 
case. However, the ratio of the two forces in the parallel and perpendicular configurations at 
small separation has a magnitude greater than 50 for R = 0.632 cm = 2 a + 0.002 cm, which 
is much larger than the factor of 2 expected from the dipole-dipole approximation. Further 
examples of this ratio are given in Table I for various separations. 

In Figs. El and I3J we test the effect of different inclusion dielectric constants, by calculating 
the force between two identical spheres, each of radius a and dielectric constant ei, in a host 
of dielectric constant Ch = 1. We plot the radial component of this force, for both the parallel 
and perpendicular configurations, as a function of ei, for two different separations between 
the spheres: R = 2a + 0.01 mm and R = 2a + 10.00 mm, where we again use a = 3.15 mm. 
In the second case, the forces are very close to the dipole-dipole predictions. In the first case, 
the forces exhibit a large departure from the predictions of the dipole-dipole interaction, and 
this departure becomes greater as ei deviates more and more from unity. 

Next, we consider an example in which the dielectric functions of both the inclusion and 
the host depend on frequency. Specifically, we choose 

. 47r(Ti 

ei = eio + i , (34) 
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and 



eh = eho + ^ 



. 47rcrh 



(35) 



where i is the imaginary unit, eio and eho are the dielectric constants of the inclusion and 
the host, and ai and ah are their conductivities, assumed frequency-independent. The time- 
averaged forces are now calculated from the generalization of eq. ()2(')|) to finite frequencies, 
namely, eq. (j^ . 

We have chosen to use parameters given by Ref. [l^, in a recent experimental study. 
These are listed in Table II. However, as discussed further below, it is possible that the 
experimentally measured forces include effects beyond the purely electrostatic interactions 
included in our model (such as spatially dependent conductivities of the host fluid). There- 
fore, our numerical results should again be considered as model calculations, not necessarily 
applicable to the specific experiments of Ref. In all cases, we assume that the two 

spherical inclusions are identical, with a dielectric constant and conductivity characteristic 
of SrTiOcj. For the host fluid, we have considered the various materials used in the mea- 

on force; we have checked this by recalculating the forces with the conductivity set equal to 
zero, and obtained the same results.) 

In Fig. El (a) and (b), we show the radial component of the calculated time-averaged 
force on a sphere of SrTiOa at R in the parallel and perpendicular geometries, for the host 
materials of silicone oil and castor oil. In both cases, we assume spheres of radius a = 
3.15 mm, intersphere spacing 6 = 0.01mm and applied electric field Eq = 25.2 V/ mm, as in 
Ref. . The magnitude of force decreases with increasing frequency, but rapidly converges 
to a constant value in both cases. The sign of the force is negative in (a), corresponding to 
an attractive force, and positive (repulsive) in (b). 

If these were strictly dipole-dipole forces, the time-averaged force on the sphere at R 
would be given by the generalization of eq. (j33|l to complex dielectric functions and eh 7^ 1, 
namely 



Thus, in particular, the magnitude of the force in the parallel case would be twice as large 
as that in the perpendicular case, as in our previous examples. However, in Fig. this force 
ratio is about 50. This difference occurs, as in Fig. |21 because of the very small separation 




(36) 
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{S = 0.01mm), which corresponds to a very short-ranged interaction. In the long range 
hmit (i? ^ a), our calculated magnitude ratio agrees well with the dipole-dipole prediction, 
as discussed further below. This short-distance deviation from dipole-dipole forces is similar 
to that seen in Figs. 

Fig. El also shows that there is a substantial difference between the forces for silicone oil 
and castor oil hosts. This difference is due almost entirely to the difference in the static 
dielectric constants of these two hosts: the effect of the finite conductivity disappears by 
about 10 Hz in both cases, whereas the difference between the forces persists to much higher 
frequencies. 

The calculated time-averaged force between spheres of SrTiOs in a silicone oil host is 
plotted versus separation in Fig. El at a frequency of 50 Hz. In order to see the effects of 
a finite host conductivity, we include this conductivity in Figs. El (a) and El (b) but not in 
El (c) or El (d). We also set the conductivity of the sphere equal to zero in (c) and (d). 
Clearly, the host conductivity has very little influence on the forces at this frequency. For 
comparison, we also show the forces as calculated in the dipole-dipole approximation. As 
can be seen, there is very little difference between the two except for R <~ 1.5 cm. Even at 
such small spacings, the deviation from the dipole-dipole force is much larger for the parallel 
than the perpendicular configuration. At a spacing of 0.01mm, the calculated ratio of force 
magnitudes in the parallel and perpendicular configurations exceeds a factor of 100. 

At sufficiently high host conductivity, our model predicts that the force between spheres 
changes sign as a function of frequency. This trend is shown in Fig. (3 for a separation 
of 5 = 0.01mm between spheres. The host materials used here are ethyl benzoate, ethyl 
salicylate, and methyl salicylate, all of which have much greater conductivities than silicone 
oil. The sign change is due mainly to the greater conductivities, not the differences in static 
dielectric constants. To check this point, we recalculated the points of Fig. [71 assuming the 
same value of the real part of the dielectric constant for all three host materials; we found that 
the time-averaged forces changed sign at the same frequencies as in Fig. |3 Mathematically, 
the origin of the sign change is, of course, the dependence of the variable s in eqs. (jHI) and 
fl29|) on the host conductivity. 

The time-averaged force for this separation ranges from about +1.5 to —1.5 dynes for 
the parallel case, depending on the frequency, and from about +0.5 to —3.0 dynes for the 
perpendicular case. At high frequencies, the force approaches —1.0 dyne for the parallel case. 
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whatever the host fluid is, and approaches a much smaller magnitude in the perpendicular 
case. The ratio of these forces differs greatly from the predictions of the dipole-dipole 
interaction, as expected for such a small separation. At very low frequencies, however, the 
force ratio appears to approach the dipole-dipole prediction. 

Fig. IHl shows the frequency dependence of the time-averaged force between two spheres of 
SrTiOs for silicone oil and N2 hosts. Both the spacings 6 between the two spheres and the 
electric fleld Eq are larger than those for Fig. El they are given in the legends of each Figure. 
We chose these values for the parameters because they are used in the measurements of 

n 

Ref . . Evidently, the force between the two spheres is stronger when the two spheres are 
immersed in a liquid host than in a gas, all the other parameters of the forces being held 
constant. This behavior can be understood even in the dipole-dipole limit: it is due to the 
dependence of the force on eh as in eq. (jHUj) . Also, the low-frequency forces in Fig. |S1 (a) 
and (b) and especially (c) and (d) depend more weakly on frequency than those in Fig. |31 
Another point is that, even though the intersphere spacing 5 has been increased to 0.10 and 
0.30 mm in these calculations, the calculated forces are still far from the dipole-dipole limit. 
Speciflcally, the ratio of the force magnitudes in the parallel and perpendicular geometries 
greatly exceeds the factor of two expected in the dipole-dipole limit. However, this ratio is 
smaller than that of Fig. presumably because the intersphere separations are larger than 
in that Figure. 



IV. DISCUSSION 



The present work permits calculation of electrical forces in ER fluids in a concise closed 
form, which permits inclusion of all multipoles and all many-body forces in a simple way. In 
our approach, the forces do not need to be calculated as numerical derivatives; instead, we 
give explicit analytical expressions for these derivatives, in terms of a pole spectrum which 
characterizes the microgeometry of the material. The explicit form for the derivatives is 
somewhat reminiscent of the Hellmann-Feynman description of quantum-mechanical forces 
in electronic structure theory, but differs from it in the important respects. 

One striking feature of the present formalism is that it allows for the calculation of 



frequency-dependent forces in a simp 



cussed in previous work 







11 
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e closed form. Although such forces have been dis- 



2J], the present approach is relatively simple and more 
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general, and places both zero and finite frequency forces within the same formalism. In our 
numerical work, we find that these forces can even change sign as a function of frequency. 
Such frequency- dependence is, of course, also present in the long-range (dipole-dipole) limit 
treated by others in the previous work, but it is even more apparent in the present study. 

Although in the present work calculations have been carried out explicitly for two-body 
interaction, they can readily be extended to three-body (or multi-body) forces. The general 
equation (pUj) or ()29|) can be used to find the force on a sphere, no matter how many 
particles are contained in the suspension. Indeed, such multi-body forces are very likely 
to play important roles in dense suspensions, where they could possibly lead to "bond- 
angle-dependent" forces analogous to angle dependent interatomic elastic forces in liquid 
and solid semiconductors. Likewise, the calculations could be readily extended to more 
complex particles (e. g., hollow spherical shells), and to non-spherical particles, provided 
that the requisite pole spectra and matrix elements can be calculated. Also, although we 
have restricted our calculations in this paper to the radial component of the interparticle 
forces, other components can be straightforwardly computed. Finally, the present formalism 
can be immediately extended to the important case of magnetorheological fluids. For such 
fluids, eq. ()26|) or (j^^ for the force would continue to be valid, provided that ei and eh are 
replaced by /ii and fi^. 

Our calculated frequency- dependent forces, obtained using parameters quoted for SrTiOs 
spheres in a conducting host, may appear to disagree with those obtained in Ref. |3l at 
close spacing. One possible explanation for this discrepancy is that the host fluid does not 
exhibit its usual bulk conductivity when two highly polarizable spheres are placed in it 
in close proximity. Instead, there could well be non-linear screening effects of the Debye- 
Hiickel type^, which would mean that the picture of a two-component composite is simply 
not appropriate in this regime. In support of this hypothesis, we note that the reported 
experimental forces are still frequency- dependent at high frequencies, while the complex 
dielectric functions of both host and sphere should be nearly frequency-independent in this 
regime, leading to a frequency- independent force in this range. 

The present method could readily be combined with standard molecular dynamics ap- 
proaches to compute dynamical properties of electrorheological (or magnetorheological) flu- 
ids. Specifically, one could carry out molecular dynamics (MD) calculations, following the 
approach of several authors^, Q, 42, 13, 13- such approaches, the force on a given 
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sphere is typically expressed as the sum of a hard-sphere repulsion, a viscous force, and 
an electrostatic force. The first two of these forces would be the same as in the previous 
MD studies, but the third would be calculated using the present method, rather than the 
dipole-dipole force generally used in most other MD studies. It would be of great interest 
to see how such quantities as viscous relaxation time would be affected by using our forces 
in these calculations. In addition to such calculations, one could study minimum-energy 
configurations of dielectric suspensions in an applied electric field, based on the forces cal- 
culated using the methods outlined here. Many such studies can already be found in the 
literature (see, e. g., Ref. 3] or It would be of interest to extend the present approach 

to calculating minimum-energy configurations including non-dipolar forces, as outlined in 
the present work. 
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TABLES 



i?(cm) 


Force ratio 


(i?-2a)/(2a) 


0.630 


602.3 


0.0000 


0.631 


88.5 


0.0016 


0.632 


52.5 


0.0032 


0.633 


38.8 


0.0048 


0.634 


31.5 


0.0063 


0.635 


26.8 


0.0079 


0.636 


23.5 


0.0095 


0.637 


21.1 


0.0111 


0.638 


19.2 


0.0127 


0.639 


17.7 


0.0143 


0.640 


16.5 


0.0159 



TABLE I. The ratios of the magnitudes of the forces between two identical spheres in the paraUel and 
perpendicular configurations, calculated at several small separations and assuming ei — 10^, eh = 1, a = 
3.15 mm, Eq = 25.2 V/mm, and oj = Q. The force is attractive in the parallel configuration, repulsive in the 
perpendicular configuration. 



Material 


Dielectric constant 


Conductivity 


SrTiOg 


249.0 


2.0 X 10-*^ 


Silicone oil 


2.54 


1.0 X 10-13 


Castor oil 


4.20 


1.0 X 10-13 


Ethyl benzoate 


5.45 


5.0 X 10-s 


Ethyl salicylate 


8.65 


1.0 X 10-^ 


Methyl salicylate 


9.46 


6.0 X 10-^ 


A^2 gas 


1.00058 






TABLE IL Parameters for the calculations shown in Figs. [SHHl The columns denote the material, the real 
part of its dielectric constant, and its conductivity (in S/m). All except for SrTiOs are used as host materials 
in the suspensions. 
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(a) 





(b) 



FIG. 1: Geometry considered in most of our calculations: Two identical spheres of radius a are located at 
the origin and sX z = R, and are contained in a host material. 5 is the surface-to-surface distance between 

the two spheres. The complex dielectric function of the spheres is and that of the host material is 

eh(w). A spatially uniform electric field is applied in the z direction in (a) and in the x direction in (b). 

FIGURES 
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FIG. 2: Magnitude of the radial component of the force at zero frequency between two identical spheres 
of radius a, with ej = 10^, = 1, plotted as a function of sphere separation, for electric field parallel to 
axis between spheres (^Eq — 0), and field perpendicular to that axis (0Eo — 7''/2). Note the logarithmic 
scale on the vertical axis. In both cases, we assume sphere radii of 3.15 mm, and an electric field of strength 
25.2 V/ mm, as in Ref. 15]. The force in the parallel field case is negative (attractive) while that in the 
perpendicular field case is positive (repulsive). In this and the following two plots, the force is calculated at 
zero frequency. 




0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 



logio Ei logio Ei 

FIG. 3: The radial component of the force at zero frequency between two identical spheres of dielectric 
constant ej, radius a = 3.15 mm, in a host of dielectric constant eh = 1, at an intersphere spacing (surface- 
to-surface separation) of 0.01 mm, plotted as a function of e^, for (a) electric field parallel to the axis between 
spheres, and (b) field perpendicular to that axis. We assume an electric field of strength 25.2 V/mm. Negative 
and positive forces denote attractive and repulsive forces, respectively. 
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FIG. 4: Same as Fig. but for an intersphere spacing S = 10.00 mm. 
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FIG. 5: The radial component of the time-averaged force between two identical spheres of SrTiOs, plotted 
as a function of frequency for the host materials of silicone oil and castor oil, respectively. For both cases 
we use 5 = 0.01 mm, a — 3.15 mm, and Eq = 25.2 V/mm. The electric field is parallel to the line connecting 
the spheres in (a) and perpendicular to that line in (b). A negative value denotes an attractive force. 



26 




0.63 1.00 



1.50 2.00 
R(cm) 



2.50 



3.00 



■/statvolt" 


1.0 n 

0.0 - 


dipole-dipole force 
numerical results - 




B 






(c) 


(dyne 


-1.0 - 






-2.0 - 






W 








bT 


-3.0 - 






o 

60 
O 


-4.0 - 


II 



0.63 1.00 



1.50 2.00 
R(cm) 



2.50 



3.00 



S 



a 



a 



— . 

ft- 




0.63 1.00 



0.63 1.00 



1.50 2.00 
R(cm) 



2.50 



3.00 



dipole-dipole force 
numerical results 




1.50 2.00 
R(cm) 



2.50 



3.00 



FIG. 6: (a) and (b) : Magnitude of the radial component of the time-averaged force between two identical 
spheres of SrTiOs, divided by Eq. Also plotted is the corresponding quantity in the dipole-dipole approxi- 
mation. Both are plotted on logarithmic scale as a function of separation R for a host material of silicone 
oil and a fixed frequency of 50 Hz. The spheres have radii 3.15 mm. The electric field is parallel to the line 
between the two spheres in (a) and perpendicular to that line in (b). (c) and (d) : Same as (a) and (b) 
except that the conductivities of the spheres and the host are set equal to zero in these calculations. The 
forces are attractive in (a) and (c), repulsive in (b) and (d). 
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(b) 
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FIG. 7: The radial component of the time-avcragcd force between two identical spheres of SrTiOs separated 
by R, plotted as a function of frequency for the host materials of ethyl benzoate, ethyl salicylate, and 
methyl salicylate, respectively. The electric field is parallel to the line connecting the two spheres in (a) and 
perpendicular to that line in (b). In all cases, 6 = 0.01mm, a = 3.15 mm, and Eq = 25.2 V/ mm. 
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FIG. 8: The radial component of the time-averaged force between two identical spheres of SrTiOs separated 
by R, plotted as a function of frequency for host materials consisting of silicone oil [(a) and (b)] and N2 [(c) 
and (d)], with gap spacings S = 0.10 mm and S = 0.30 mm. The applied electric field is £^0 = 71.3 V/ mm 
and a = 3.15 mm for all the cases. The electric field is parallel to the line between two spheres in (a) and 
(c) and perpendicular to that line in (b) and (d). 
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